### This project replicates the analyses reported in Toshkov and Krouwel (2022) 'Beyond the U-Curve: Citizen Preferences on European Integration in Multidimensional Political Space' in European Union Politics 
### Part III: This script prepares and analyzes the data from the 2019 EP elections in Italy
# Clean up VAA data and make variables -------------------------------------------------------------------------
## This part of the script shows the preparation of the datafile that is shared from the original datafile that is proprietary
# Load the data
it.all <- read_sav('./data in/it ep 2019/Italy EU2019 AN.sav')

# Recode variables and filter
it.all <- it.all %>%
  mutate(lr = popup_answer_1059 - 2,
         proeu = (4-popup_answer_1061) - 2,
         progr = popup_answer_1060 - 2,
         ep2019.vote.intent = recode (popup_answer_1058, 
                                      '0' = 'PD', '1' = 'M5S', 
                                      '3' = 'Lega', '2' = 'Forza Italia', 
                                      '5'= 'La Sinistra', '4'= 'PPI',
                                      '6' = 'Fratelli d`Italia', '7'='Europa Verde',
                                      '8' = '+ Europa', 
                                      '12' = 'Partito Comunista', '13' = 'Partito Pirata',
                                      .default = 'Other'),
        tk2018.vote = recode (popup_answer_1056, 
                               '0' = 'M5S', '1' = 'Lega', 
                               '2' = 'PD', '3' = 'Forza Italia', 
                               '4' = 'Fratelli d`Italia', '5' = 'Liberi e Uguali',
                               '6' = '+ Europa', '8' = 'Potere al Popolo!',
                               '17' = 'No vote', .default = 'Other'),
         themes = case_when ((relevant_themes_182 == 1 ~ 'immigration'),  
                             (relevant_themes_183 == 1 ~ 'eu'),
                             (relevant_themes_185 == 1 ~ 'security'), 
                             (relevant_themes_186 == 1 ~ 'economic policy'),
                             (relevant_themes_187 == 1 ~ 'redistribution'), 
                             (relevant_themes_188 == 1 ~ 'moral values')),
         type = 'voter',
         eu.bad = statement_answer_5031 - 3,
         euro.out = statement_answer_5032 - 3,
         eu.fp = statement_answer_5034 - 3,
         eu.elect.com = statement_answer_5039 - 3,
         eu.army = statement_answer_5040 - 3,
         tax.multinats = statement_answer_5041 -3,
         eu.too.far = statement_answer_5042 -3,
         eu.asylum.solidarity= statement_answer_5043 -3,
         trade.china = statement_answer_5044 -3,
         eu.citizens= statement_answer_5033 -3,
         islam.bad = statement_answer_5049 -3,
         free.abortion= statement_answer_5050 -3,
         immigrants.adapt = statement_answer_5052 -3,
         defend.property = statement_answer_5054 -3,
         indep.judges = statement_answer_5055 -3,
         more.citizenship = statement_answer_5056 -3,
         no.support.homo.couples = statement_answer_5057 -3,
         redistribute = statement_answer_5059 -3,
         state.econ.out = statement_answer_5060 -3,
         fire.easy = statement_answer_5062 -3,
         more.market.health = statement_answer_5063 -3,
         protect.environment = statement_answer_5064 -3,
         austerity = statement_answer_5065 -3,
         eu.allow.nat.aid = statement_answer_5067 -3,
         support.poor = statement_answer_5069 -3,
         reduce.public.debt = statement_answer_5070 -3,
         flat.tax = statement_answer_5071 -3,
         homo.rights = statement_answer_5046 -3,
         weed.legal = statement_answer_5047 -3,
         limit.privacy.fight.crime = statement_answer_5048 -3,
         ) %>%
  filter (country == 'Italy', amount_visits == 1,
          seconds_on_statements > 30*3, seconds_on_statements < 30*30) %>%
  select (lr:ncol(.))

it.all$lr <- ifelse (it.all$lr>=3, NA, it.all$lr)
it.all$proeu <- ifelse (it.all$proeu<=-3, NA, it.all$proeu)
it.all$progr <- ifelse (it.all$progr>=3, NA, it.all$progr)

# remove those who answered everything with Neutral
it.all$zeroes <- it.all %>% 
  select (eu.bad:limit.privacy.fight.crime)  %>%
  abs () %>%
  rowSums(na.rm=T)

it.all <- it.all %>% filter (zeroes!=0) %>% select (-zeroes)



it.p <- it.p [-c(5,10,14,15,16),]

it.all <- full_join(it, it.p)

# create deductively derived positions on scales
it.all <- it.all %>%
  mutate (obj.lr = (- tax.multinats - redistribute + state.econ.out + fire.easy + more.market.health + austerity - eu.allow.nat.aid + reduce.public.debt + flat.tax - support.poor)/10,
          obj.cp = (- islam.bad	+ free.abortion + indep.judges + more.citizenship - no.support.homo.couples - immigrants.adapt - defend.property - protect.environment	+
                      homo.rights	+ weed.legal - limit.privacy.fight.crime)/11,
          obj.eu = (- eu.bad - euro.out + eu.fp + eu.elect.com + eu.army - eu.too.far - eu.citizens)/7
          )

it.all$id = seq.int(nrow(it.all)) #index for merging later

# save
save (it.all, file = './data out/it2019_replication.RData')

# Load and subset the data (start replication here) -----------------------------------------------------------
load (file = './data out/it2019_replication.RData')

# PCA and Factor analysis ----------------------------------------------------------------
it.s <- it.all %>% select(id, eu.bad:limit.privacy.fight.crime) %>%
  filter(complete.cases(.))

# get polychoric covariance
tk.poly<-polychoric(select(it.s, -id))

# determine number of factors
ev <- eigen(cor(select(it.s, -id)))
ev.poly <- eigen(tk.poly$rho) #with the polychronic correlations
# number of factors

ap <- parallel(subject=nrow(it.s),var=ncol(select(it.s, -id)),rep=10,quantile=.95)
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea)
plotnScree(nS)

# principal components (Table A12)
pca.1 <- principal(select(it.s, -id), nfactors=5, rotate='none', cor='poly', scores=TRUE)
print(pca.1, digits=2, cut=.33, sort=T)

pca.out<-as.data.frame(unclass(pca.1$loadings))
pca.out <- pca.out %>%
  mutate_at(vars(starts_with('PC')), ~replace(., abs(.) < 0.33, 0)) %>%
  mutate_at(vars(starts_with('PC')), funs(round(., 2))) %>%
  arrange (-abs(PC1),-abs(PC2), -abs(PC3), -abs(PC4), -abs(PC5)) %>%
  mutate_at(vars(starts_with('PC')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
pca.out <- rbind (pca.out, c('Proportion Variance', paste0(round(prop.table(pca.1$values)[1:5]*100,0), '%')))
pca.out

pca.it.table = htmlTable::htmlTable(pca.out, rnames = FALSE,
                                    col.columns = c('white','lightgrey'),
                                    tfoot = 'Note: Loadings smaller than +/-0.41 are not printed',
                                    align = c('l','r','r','r','r','r'),
                                    align.header = c('l','r','r','r','r','r'),
                                    header = c('Policy issue','PC1','PC2','PC3','PC4','PC5'),
                                    caption = "Table X. PCA results from Italy, 2019")
print(pca.it.table, type = "html", file = './tables/pca_it_table.html')

# factor analysis with polychoric correlations and oblique rotation (Table 3)
fa.1 <- fa (select(it.s, -id), nfactors=5, rotate='Promax', cor='poly', scores=TRUE)
print(fa.1, digits=2, cut=.27, sort=T)

fa.out<-as.data.frame(unclass(fa.1$loadings))
fa.out <- fa.out %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.27, 0)) %>%
  mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>%
  arrange (-abs(MR1),-abs(MR4), -abs(MR3), -abs(MR2), -abs(MR5)) %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 
fa.out <- rbind (fa.out, c('Proportion Variance', paste0(round(fa.1$Vaccounted[2,]*100,0), '%')))

fa.out

fa.it.table = htmlTable::htmlTable(fa.out, rnames = FALSE,
                                   col.columns = c('white','lightgrey'),
                                   tfoot = 'Note: Loadings smaller than +/-0.31 are not printed',
                                   align = c('l','r','r','r','r','r','r'),
                                   align.header = c('l','r','r','r','r','r','r'),
                                   header = colnames(fa.out),
                                   caption = "Table X. Factor Analysis results from Italy, 2019")
print(fa.it.table, type = "html", file = './tables/fa_it_table.html')

# Run the same factor model but only on the subset with party affiliation (Table A13)
it.s2 <- it.all %>%
  select (id, lr, progr, proeu, obj.lr, obj.cp, obj.eu, eu.bad:limit.privacy.fight.crime) %>%
  na.omit 

# factors with fa.poly and promax rotation
fa.2 <- fa (select(it.s2, eu.bad:limit.privacy.fight.crime), nfactors=5, rotate='Promax', cor='poly', scores=TRUE)
print(fa.2, digits=2, cut=.33, sort=T)

fa2.out<-as.data.frame(unclass(fa.2$loadings))
fa2.out <- fa2.out %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.27, 0)) %>%
  mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>%
  arrange (-abs(MR1),-abs(MR3), -abs(MR4), -abs(MR2), -abs(MR5)) %>%
  mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>%
  mutate ( issue = rownames(.)) %>%
  relocate (issue) 

fa2.out <- rbind (fa2.out, c('Proportion Variance', paste0(round(fa.2$Vaccounted[2,]*100,0), '%')))

fa2.it.table = htmlTable::htmlTable(fa2.out, rnames = FALSE,
                                    col.columns = c('white','lightgrey'),
                                    tfoot = 'Note: Loadings smaller than +/-0.31 are not printed',
                                    align = c('l','r','r','r','r','r','r'),
                                    align.header = c('l','r','r','r','r','r','r'),
                                    header = colnames(fa2.out),
                                    caption = "Table SS6. Factor Analysis results from Italy, 2019")
print(fa2.it.table, type = "html", file = './tables/fa2_it_table.html')

# assign scores to cases
it.s <- bind_cols(it.s, data.frame(pca.1$scores), data.frame(fa.1$scores))

it.all <- left_join (it.all, select(it.s, id, PC1:MR5), by='id')

# Subset -------------------------------------------------
it.s <- it.all %>% 
  filter (type=='voter') %>% 
  select ('lr','progr','proeu') %>% 
  filter(complete.cases(.))

# Self-placement analysis: 3D plots -------------------------------------------------
#prepare data
elevation.df = data.frame(x = it.s$lr, y = it.s$progr, z = it.s$proeu)
elevation.loess = loess(z ~ x*y, data = elevation.df, degree=2, span=0.8)
xmin<- -2
xmax<- 2
ymin<- -2
ymax<- 2

elevation.fit = expand.grid(list(x = seq(xmin, xmax, 0.1), y = seq(ymin, ymax, 0.1)))
z = predict(elevation.loess, newdata = elevation.fit)
elevation.fit$Height = as.numeric(z)

#prepare matrix with x and y and rows/cols and z as the cell values
u1<-unique(elevation.fit$x)
u2<-unique(elevation.fit$y)
temp.data<-matrix(NA, nrow=length(u1), ncol=length(u2))

for (i in 1:length(u1)){
  for (j in 1:length(u2)){
    temp.data[i,j]<-elevation.fit$Height[elevation.fit$x==u1[i] & elevation.fit$y==u2[j]]
  }
}

ix <- seq(1, nrow(temp.data), length.out = 11)
iy <- seq(1, ncol(temp.data), length.out = 11)

# 3D plots in black and white (Figure 6, left panel)
png('./figures/it 2019 eu3.png', width=900*2, height=700*2, res=96)
par(mfrow = c(1, 1), mar = c(0, 0, 0, 0))
persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, 
        shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")
dev.off()

# 4 3d plots together (Figure A8, top)
png('./figures/it 2019 eu6.png', width=900*2, height=700*2, res=96)
par(mfrow = c(2, 2), mar = c(3, 3, 3, 3))

persp3D(z = temp.data, theta = 55, phi = 20, facets = T, curtain = TRUE, 
        shade=0.8, colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU", box=T, bty='b2',ticktype = "detailed")

persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, 
        shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

ribbon3D(z = temp.data[ix, ], along = "y", cex.axis = 0.8, 
         curtain = T, space = 0.7, theta = 25, phi = 5, shade=0.2, colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU", 
         box=T, bty='b2', ticktype = "detailed")

hist3D(z = temp.data[ix,iy], theta = 55, phi = 20, shade=0.2, colkey = FALSE, box=T, bty='b2',ticktype = "detailed",
       xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

dev.off()

# now with the rayshader package (Figure A8, bottom)
temp.data2 = temp.data*300

raymat = ray_shade(temp.data2)
ambmat = ambient_shade(temp.data2)

png('./figures/raysahder_it.png', width=1000, height=800, res=96)

temp.data2 %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(temp.data2), color="desert") %>%
  add_shadow(ray_shade(temp.data2,zscale=3,maxsearch = 300),0.5) %>%
  add_shadow(ambmat,0.5) %>%
  plot_3d(temp.data2,zscale=35,fov=0,theta=-197,zoom=0.75,phi=10, windowsize = c(1000,800),
          water=TRUE, waterdepth = 0, wateralpha = 0.5,watercolor = "lightblue",
          waterlinecolor = "white",waterlinealpha = 0.5)

render_label(temp.data2, x=1,y=1, z=510, zscale=35,relativez=FALSE, text = "Left, Conservative",textsize = 1.7,linewidth = 2,freetype = FALSE,
             linecolor = 'darkred', textcolor = "darkred") 
render_label(temp.data2, x=41,y=1, z=510, zscale=35,relativez=FALSE,
             linecolor = 'darkblue', textcolor = "darkblue", text = "Right, Conservative",textsize = 1.7,linewidth = 2,freetype = FALSE)
render_label(temp.data2, x=1,y=41, z=510, zscale=35,relativez=FALSE,
             linecolor = 'darkred', textcolor = "darkred", text = "Left, Progressive",textsize = 1.7,linewidth = 2,freetype = FALSE)
render_label(temp.data2, x=41,y=41, z=510, zscale=35, relativez=FALSE,
             linecolor = 'darkblue',textcolor = "darkblue", text = "Right, Progresive",textsize = 1.7,linewidth = 2,freetype = FALSE)

#render_label(temp.data2, x=26,y=101, z=1400, zscale=35, relativez=FALSE, dashed= TRUE,
#             linecolor = 'black', textcolor = "black", text = "Peak of EU support",textsize = 1.5,linewidth = 2,freetype = FALSE, offset=0)
#render_label(temp.data2, x=85,y=15, z=300, zscale=35, relativez=FALSE, dashed= TRUE,
#             linecolor = 'black', textcolor = "black", text = "Lake of Eurosceptics",textsize = 1.5, linewidth = 2,freetype = FALSE)

render_snapshot()
dev.off()

# Additional 2D and 3D density plots --------------------------------------------------

# 3D histogram with the density mapped as color
ix <- seq(1, nrow(temp.data), length.out = 6)
iy <- seq(1, ncol(temp.data), length.out = 6)
smooth <- MASS::kde2d(it.s$lr, it.s$progr, lims = c(-2, 2, -2, 2), n = 6)$z 

# Figure A9
png('./figures/it 2019 dens eu.png', width=900*2, height=700*2, res=96)
par(mfrow = c(2, 2), mar = c(3, 3, 3, 3))
color = rev(rainbow(10, start = 0/6, end = 4/6))

hist3D(z = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, colvar = smooth, col=color, colkey = FALSE,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

hist3D(z = temp.data[ix,iy], theta = 35 + 180, phi = 20, shade=0.2, colvar = smooth, col=color,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

hist3D(colvar = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, z = smooth, col=color, colkey = FALSE,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density")

hist3D(colvar = temp.data[ix,iy], theta = 35 + 180, phi = 20, shade=0.2, z = smooth, col=color,
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density")
dev.off()

# Subjective placement and objective positions: correlation plots (Figure 4) ----------------------------------
png('./figures/it 2019 corplot2.png', width=900, height=900, res=96)

ggpairs(
  it.all [it.all$type=='voter', c('obj.lr', 'obj.cp', 'obj.eu', 'lr','progr','proeu')],
  upper = list(continuous = "cor"),
  #upper = list(continuous = "density"),
  lower = list(continuous = "trends"),
  diag  = list(continuous = "barDiag"),
  axisLabels = 'show',
  columnLabels = c('Obj. Left-Right','Obj. Cons.-Progr.','Obj. Anti-Pro EU',
                   'Left-Right','Conservative-Progressive','Anti-Pro EU')) + 
  theme_bw()
dev.off()

# Objective positions: 3D plots ----------------------------------
it.s2 <- it.all %>% 
  filter (type=='voter') %>% 
  select ('obj.lr','obj.cp','obj.eu') %>% 
  filter(complete.cases(.))

#prepare data
elevation.df = data.frame(x = it.s2$obj.lr, y = it.s2$obj.cp, z = it.s2$obj.eu)
elevation.loess = loess(z ~ x*y, data = elevation.df)

xmin<- -2
xmax<- 2
ymin<- -2
ymax<- 2

elevation.fit = expand.grid(list(x = seq(xmin, xmax, 0.1), y = seq(ymin, ymax, 0.1)))
z = predict(elevation.loess, newdata = elevation.fit)
elevation.fit$Height = as.numeric(z)

#prepare matrix with x and y and rows/cols and z as the cell values
u1<-unique(elevation.fit$x)
u2<-unique(elevation.fit$y)
temp.data<-matrix(NA, nrow=length(u1), ncol=length(u2))

for (i in 1:length(u1)){
  for (j in 1:length(u2)){
    temp.data[i,j]<-elevation.fit$Height[elevation.fit$x==u1[i] & elevation.fit$y==u2[j]]
  }
}

ix <- seq(1, nrow(temp.data), length.out = 11)
iy <- seq(1, ncol(temp.data), length.out = 11)

smooth <- MASS::kde2d(it.s2$obj.lr, it.s2$obj.cp, lims = c(-2, 2, -2, 2), n = 11)$z 

# 3D plot in black and white (Figure 6, right panel)
png('./figures/it 2019 obj eu3.png', width=900*2, height=700*2, res=96)
par(mfrow = c(1, 1), mar = c(0, 0, 0, 0))
persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, 
        shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Obj. Left-Right", ylab="Obj. Conservative-Progressive", zlab="Obj. Anti-Pro EU")
dev.off()


# additional 3D plots (not included)
png('./figures/it 2019 dens eu3.png', width=900, height=700, res=96)

par(mfrow = c(2, 2), mar = c(1, 1, 1, 1))
color = rev(rainbow(10, start = 0/6, end = 4/6))

hist3D(z = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, colvar = smooth, col=color, colkey = FALSE, clab='Density',
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

hist3D(z = temp.data[ix,iy], theta = 35 - 90, phi = 20, shade=0.2, colvar = smooth, col=color, clab='Density',
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU")

hist3D(colvar = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, z = smooth, col=color, colkey = FALSE, clab='Anti-Pro EU',
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density")

hist3D(colvar = temp.data[ix,iy], theta = 35 - 90, phi = 20, shade=0.2, z = smooth, col=color, clab='Anti-Pro EU',
       box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density")
dev.off()

### THE END


